library(scater)
library(scran)
library(TSCAN)
library(slingshot)
library(tradeSeq)
library(Seurat)
library(tidySingleCellExperiment)
library(tidyverse)

# example data --------
sce_10x <- read_rds('CRC-I/data/crc10x_singler_ref.rds')

sce_10x |> count(Global_Cluster)

sce_10x_myl <- sce_10x |> filter(str_detect(Global_Cluster, 'Myeloid'))

sce_10x_myl <- runPCA(sce_10x_myl, ntop = 2000,
                  BSPARAM = BiocSingular::RandomParam())

colLabels(sce_10x_myl) <- sce_10x_myl$Sub_ClusterID

plotReducedDim(sce_10x_myl, dimred = 'PCA', colour_by = 'Sub_ClusterID')

sce_10x_myl <- runUMAP(sce_10x_myl)

# quickPseudotime to get MST -----------
pseudo.all <- sce_10x_myl |> TSCAN::quickPseudotime(use.dimred = 'PCA')

pseudo.all |> glimpse()

pseudo.all$connected$UMAP |>
  ggplot(aes(UMAP1, UMAP2, group = edge)) +
  geom_path()

## get pseudotime -----------
sce_10x_myl <- pseudo.all$ordering |>
  averagePseudotime() |>
  as_tibble(rownames = '.cell') |>
  rename('ptime' = value) |>
  left_join(sce_10x_myl, y = _)

sce_10x_myl |>
  plotUMAP(colour_by='ptime', text_by="label", text_colour="red") +
  geom_line(data=pseudo.all$connected$UMAP, aes(UMAP1, UMAP2, group=edge))

## introduce outgroup to identify alt MST --------
pseudo.og <- sce_10x_myl |>
  TSCAN::quickPseudotime(use.dimred = 'PCA', outgroup = TRUE)

# primary plot of igraph
plot(pseudo.og$mst)

## MST via MNN instead of cluster centroids -----------
pseudo.mnn <- sce_10x_myl |>
  TSCAN::quickPseudotime(use.dimred = 'PCA', dist.method = 'mnn', outgroup = TRUE)

sce_10x_myl <- pseudo.mnn$ordering |>
  averagePseudotime() |>
  as_tibble(rownames = '.cell') |>
  rename('ptime.mnn' = value) |>
  left_join(sce_10x_myl, y = _)

sce_10x_myl |>
  plotUMAP(colour_by='ptime.mnn', text_by="label", text_colour="red") +
  geom_line(data=pseudo.mnn$connected$UMAP, aes(UMAP1, UMAP2, group=edge))

## TSCAN MST in slingshot method ------
pseudo.ss <- sce_10x_myl |>
  TSCAN::quickPseudotime(use.dimred = 'PCA', dist.method = 'slingshot')

sce_10x_myl <- pseudo.ss$ordering |>
  averagePseudotime() |>
  as_tibble(rownames = '.cell') |>
  rename('ptime.ss' = value) |>
  left_join(sce_10x_myl, y = _)

sce_10x_myl |>
  plotUMAP(colour_by='ptime.ss', text_by="label", text_colour="red") +
  geom_line(data=pseudo.ss$connected$UMAP, aes(UMAP1, UMAP2, group=edge))

# principal curves by native slingshot -----
# slingshot with small cluster cause error!
major_myl <- sce_10x_myl |>
  count(Sub_Cluster) |>
  slice_max(n, n = 6) |>
  pull(Sub_Cluster)

sce_10x_myl6 <-
  sce_10x_myl |>
  filter(Sub_Cluster %in% major_myl)

## slingshot on clusters --------
sce_myl_ss <- sce_10x_myl6 |>
  slingshot(reducedDim = 'PCA',
            cluster = colLabels(sce_10x_myl6))

ss_ptime <- sce_myl_ss |>
  embedCurves("UMAP") |>
  slingCurves(as.df = T) |>
  arrange(Order)

ss_ptime |> head()

sce_myl_ss |>
  plotUMAP(colour_by="slingPseudotime_1") +
  geom_path(data=ss_ptime, aes(UMAP1, UMAP2, group = Lineage))

## speed large data by set approx_points -------
sce_myl_ss_approx <- sce_10x_myl6 |>
  slingshot(reducedDim = 'PCA',
            cluster = colLabels(sce_10x_myl6),
            approx_points = 100)

## add an omega outgroup ---------
sce_myl_ss_omega <- sce_10x_myl6 |>
  slingshot(reducedDim = 'PCA',
            cluster = colLabels(sce_10x_myl6),
            omega = TRUE)

ss_ptime <- sce_myl_ss_omega |>
  embedCurves("UMAP") |>
  slingCurves(as.df = T) |>
  arrange(Order)

sce_myl_ss_omega |>
  plotUMAP(colour_by="Sub_Cluster") +
  geom_path(data=ss_ptime, aes(UMAP1, UMAP2, group = Lineage))

# find marker genes along pseudotime ----------
marker.pseudo <- sce_10x_myl |>
  TSCAN::testPseudotime(sce_10x_myl$ptime)

marker.pseudo <- marker.pseudo |>
  as_tibble(rownames = 'gene') |>
  filter(FDR < .05)

## upregulated along pseudotime
marker.pseudo.up <- marker.pseudo |>
  slice_max(logFC, n = 10) |>
  pull(gene)

sce_10x_myl |>
  plotExpression(features = marker.pseudo.up,
                 x = 'ptime',
                 colour_by = 'label')

## downregulated along pseudotime
marker.pseudo.down <- marker.pseudo |>
  slice_min(logFC, n = 10) |>
  pull(gene)

sce_10x_myl |>
  plotExpression(features = marker.pseudo.down,
                 x = 'ptime',
                 colour_by = 'label')

## more genes on heatmap
marker.pseudo.up50 <- marker.pseudo |>
  slice_max(logFC, n = 50) |>
  pull(gene)

sce_10x_myl |>
  plotHeatmap(features = marker.pseudo.up50,
                 order_columns_by = 'ptime',
                 colour_columns_by = 'label',
                 center = TRUE)

## use GAM instead of lm to identify DEG ------
sce_10x_myl_gam <- sce_myl_ss_omega |>
  fitGAM(genes = marker.pseudo.up50)

gam_res <- sce_10x_myl_gam |>
  patternTest()

gam_res |>
  as_tibble(rownames = 'gene') |>
  filter(!is.na(pvalue))

# find root by entropy --------
sce_10x_myl$entropy <- perCellEntropy(sce_10x_myl)

# generally lower entropy indicates higher differentiation
sce_10x_myl |>
  ggplot(aes(Sub_Cluster, entropy)) +
  geom_violin() +
  coord_flip() +
  stat_summary(fun = median, geom = 'point')
